home *** CD-ROM | disk | FTP | other *** search
/ HPAVC / HPAVC CD-ROM.iso / pc / FAQSYS18.ZIP / FAQS.DAT / CONREC.TXT < prev    next >
Encoding:
Text File  |  1991-05-27  |  23.9 KB  |  634 lines
  1. Text and code for a contouring algorithm, wrote it a long time ago but it is
  2. still useful. If you want the image see the June or July 1987 Byte article.
  3.  
  4. CONREC - A Contouring Subroutine
  5. Written by Paul D. Bourke (MSc.)
  6.  
  7. Introduction
  8. This article introduces a straightforward method of contouring some surface. 
  9. The  algorithm is implemented in a subroutine written in BASIC.
  10. Contouring aids in visualizing three dimensional surfaces on a two 
  11. dimensional medium (on  paper or in this case a computer graphics screen). 
  12. Two most common applications are  displaying topological features of an area 
  13. on a map or the air pressure on a weather map.  In all cases some parameter 
  14. is plotted as a function of two variables, the longitude and  latitude or x 
  15. and y axis. One problem with computer contouring is the process is usually  
  16. CPU intensive and the algorithms often using advanced mathematical 
  17. techniques.
  18.  
  19. CONREC
  20. To do contouring in software you need to describe the data surface and the 
  21. contour levels  you want to have drawn. The software given this information 
  22. must call the algorithm that  calculate sthe line segments that make up a 
  23. contour curve and then plot these line  segments on whatever graphics device 
  24. is available. 
  25. CONREC satisfies the above description, it is relatively simple to 
  26. implement, very  reliable, and does not require sophisticated programming 
  27. techniques or a high level of  mathematics to understand how it works. 
  28. The input parameters to the CONREC subroutine are as follows :
  29. - The number of horizontal and vertical data points designated iub and jub.
  30. - The number of contouring levels, nc.
  31. - A one dimensional array z(0:nc-1) that seves as a list of the contour 
  32. levels in  increasing order. (The order of  course can be relaxed if the 
  33. program will sort the  levels)
  34. - A two dimensional array d(0:iub,0:jub) that contains the description of 
  35. the data array  to be contoured. Each element of the array is a sample of 
  36. the surface being studied at a  point (x,y)
  37. - Two, one dimensional arrays x(0:iub) and y(0:jub) which contain the 
  38. horizontal and  vertical coordinates of each sample point. This allows for a 
  39. rectangular grid of samples. Diagram 0 illustates some of the above input 
  40. parameters.
  41.  
  42. The contouring subroutine CONREC does not assume anything about the device 
  43. that will be  used to plot the contours. It instead expects a user written 
  44. subroutine called VECOUT.  CONREC calls VECOUT with the horizontal and 
  45. vertical coordinates of the start and end  coordinates of a line segment 
  46. along with the contour level for that line segment. In the  simplest case 
  47. this is very similar the the usual LINE (x1,y1)-(x2,y2) command in BASIC.  
  48. See examples.
  49.  
  50. Algorithm
  51. As already mentioned the samples of the three dimensional surface are stored 
  52. in a two  dimensional real array. This rectangular grid is considered four 
  53. points at a time, namely  the rectangle d(i,j), d(i+1,j), d(i,j+1), and 
  54. d(i+1,j+1). The centre of each rectangle is  assigned a value corresponding 
  55. to the average values of each of the four vertices. Each  rectangle is in 
  56. turn divided into four triangular regions by cutting along the diagonals.  
  57. Each of these triangular planes may be bisected by horizontal contour plane. 
  58. The  intersection of these two planes is a straight line segment, part of 
  59. the the contour  curve at that contour height.
  60. Depending on the value of a contour level with respect to the height at the 
  61. vertices of a  triangle, certain types of contour lines are drawn. The 10 
  62. possible cases which may occur  are summarised below (also see table 1)
  63. a) All the vertices lie below the contour level.
  64. b) Two vertices lie below and one on the contour level.
  65. c) Two vertices lie below and one above the contour level.
  66. d) One vertex lies below and two on the contour level.
  67. e) One vertex lies below, one on and one above the contour level.
  68. f) One vertex lies below and two above  the contour level.
  69. g) Three vertices lie on the contour level.
  70. h) Two vertices lie on and one above the contour level.
  71. i) One vertex lies on and two above the contour level.
  72. j) All the vertices lie above the contour level.
  73.  
  74. In cases a, b, i and j the two planes do not intersect, ie: no line need be 
  75. drawn. For  cases d and h the two planes intersect along an edge of the 
  76. triangle and therefore line  is drawn between the two vertices that lie on 
  77. the contour level. Case e requires that a  line be drawn from the vertex on 
  78. the contour level to a point on the opposite edge. This  point is determined 
  79. by the intersection of the contour level with the straight line  between the 
  80. other two vertices. Cases c and f are the most common situations where the  
  81. line is drawn from one edge to another edge of the triangle. The last 
  82. possibility or case  g above has no really satisfactory solution and 
  83. fortunately will occur rarely with real  arithmetic. Diagram 2 summarises 
  84. the possible line orientations.
  85.  
  86. Example
  87. As a simple example consider one triangle with vertices labelled m1,m2 and 
  88. m3 with  heights 0,2 and 3 respectively (See diagram 3) To calculate where a 
  89. contour line at a  height of 1 should be drawn, it can be seen that this is 
  90. case f described above. Level 1  intersects line segment m1-m2 half the way 
  91. along and it intersects line segment m1-m3 one  third of the way along. A 
  92. line segment is drawn between these two points.
  93.  
  94. Subroutine
  95. In summary, CONREC takes each rectangle of adjacent data points and splits 
  96. it into 4  triangles after chooseing the height at the centre of the 
  97. rectangle. For each of the  triangles the line segment resulting from the 
  98. intersection with each contour plane. A  routine is then called with the 
  99. starting and stopping coordinates of the line segment. Listing 1 shows the 
  100. CONREC subroutine written in MicroSoft Basic (version 2 binary) for  the 
  101. Macintosh. An attempt is made at optimization by checking first to see if 
  102. there are  any contour levels within the present rectangle and second that 
  103. there are some contour  levels within the present triangle. The indices i 
  104. and j are used to step through each  rectangle in turn, k refers to each 
  105. contour level and m to the four triangles in each  rectangle. Diagram 1 
  106. gives some of the notation used for identifying the rectangles and  
  107. triangles in the subroutine.
  108. Note that for large arrays the whole data array need be stored in memory . 
  109. Since the  algorithm is a local one only requiring 4 points at a time, the 
  110. data for each rectangle  could be read from disk when required.
  111.  
  112. Improvements
  113. You can add features such as labeling the contour and axes, drawing the 
  114. contour lines  with different colours or in the case of a black and white 
  115. display, drawing the contour  lines with different mark-space ratio's. The 
  116. contours may be labeled with the contour  values either by software or by 
  117. adding them later by hand.
  118.  
  119.  
  120. Examples of CONREC
  121.  
  122. Example 1 is a program that uses CONREC to draw the contour of a 
  123. mathematical function.  The function is given by
  124.                                                      1 
  125.      f ( x , y ) =  sin ( ( x2 + y2 )1/2 ) +  ----------------
  126.                                              ( ( x - c )2 + y2 )1/2 
  127.  
  128. Example 2 is taken from physics. The potential a distance r from the origin 
  129. from two  charges q1 and q2, located at plus and minus c is given by
  130.                        q1      q2
  131. V( x , y ) = const * ( ---  -  --- )
  132.                        r1       r2
  133. where r1 and r2 are the distances from the point ( x , y ) to the charge q1 
  134. and q2  respectively. Contouring this function gives rise to equipotential 
  135. surfaces, Similar  functions also arise in the case of gravitational 
  136. potentials. 
  137.  
  138. Example 3 is a contour diagram of the function
  139.                                     1
  140. f( x , y ) = -------------------------------------------------
  141.            ((x2+(y-0.842)(y+0.842))2 + (x(y-0.842) + x(y-0.842))2
  142.  
  143. Example 4 is a contour diagram of cowboy hat shaped function
  144. f( x , y ) = sin ( x2 + y2 )
  145. The examples given here are written in Basic and are therefore rather slow. 
  146. Example 1  (31x31 data array and requesting 10 contour levals) took about 5 
  147. minutes to draw.
  148.  
  149. ****** Demo FORTRAN-77 program *******
  150.  
  151.       program demo
  152.       implicit none
  153. c
  154. c     The follwoing is a simplistic application of the CONREC routine.
  155. c     A mathematical function is evaluated over a regular grid of points
  156. c     on a computer raster graphics screen.
  157. c
  158. c     Paul D. Bourke
  159. c
  160.       integer*4 pxmin,pymin,pxmax,pymax
  161.       parameter (pxmin = 10,pymin = 10,pxmax = 460,pymax = 300)
  162.       integer*4 nx,ny
  163.       parameter (nx = 100,ny = 50)
  164.       integer*4 nc
  165.       parameter (nc = 10)
  166. c
  167.       real*4 d(0:nx,0:ny),x(0:nx),y(0:ny),z(1:nc)
  168.       real*4 x1,y1,x2,y2
  169.       real*4 zmax,zmin
  170.       integer*4 i,j
  171. c
  172. c     Create an artificial data surface and calculate the
  173. c     surface bounds for choosing the contour levels.
  174. c
  175.       zmin =  1.0e30
  176.       zmax = -1.0e30
  177.       do 200 i=0,nx
  178.          do 100 j=0,ny
  179.             d(i,j) = i * j
  180.             zmin = min(zmin,d(i,j))
  181.             zmax = max(zmax,d(i,j))
  182. 100      continue
  183. 200   continue
  184. c
  185. c     Set coordinates in Y array suitable for 
  186. c     automatic plotting on the graphics screen
  187. c
  188.       do 300 j=0,ny
  189.          y(j) = pymax - j * (pymax - pymin) / float(ny)
  190. 300   continue
  191. c
  192. c     Set coordinates in X array suitable for
  193. c     automatic plotting on the graphics screen
  194. c
  195.       do 400 i=0,nx
  196.          x(i) = i * (pxmax - pxmin) / float(nx) + pxmin
  197. 400   continue
  198. c
  199. c     Set a full contingent of contour levels
  200. c
  201.       do 500 i=1,nc
  202.          z(i) = i * (zmax - zmin) / (nc + 1)
  203. 500   continue
  204. c
  205. c     Draw a border around the contour plot
  206. c
  207.       x1 = pxmin
  208.       y1 = pymin
  209.       x2 = pxmax
  210.       y2 = pymax
  211.       call vecout(x1,y1,x1,y2,0.0)
  212.       call vecout(x1,y2,x2,y2,0.0)
  213.       call vecout(x2,y2,x2,y1,0.0)
  214.       call vecout(x2,y1,x1,y1,0.0)
  215. c
  216. c     Call the contouring routine
  217. c
  218.       call conrec(d,0,nx,0,ny,x,y,nc,z)
  219.       pause
  220. c
  221.       end
  222.  
  223. ****** FORTRAN-77 CODE *******
  224.  
  225. c
  226. c======================================================================
  227. c
  228. c     CONREC is a contouring subroutine for rectangularily spaced data.
  229. c
  230. c     It emits calls to a line drawing subroutine supplied by the user
  231. c     which draws a contour map corresponding to real*4data on a randomly
  232. c     spaced rectangular grid. The coordinates emitted are in the same
  233. c     units given in the x() and y() arrays.
  234. c
  235. c     Any number of contour levels may be specified but they must be 
  236. c     in order of increasing value.
  237. c
  238. c     subroutine conrec(d,ilb,iub,jlb,jub,x,y,nc,z)
  239. c     real*4 d(ilb:iub,jlb:jub)  ! matrix of data to contour
  240. c     integer ilb,iub,jlb,jub    ! index bounds of data matrix
  241. c     real*4 x(ilb:iub)          ! data matrix column coordinates
  242. c     real*4 y(jlb,jub)          ! data matrix row coordinates
  243. c     integer nc                 ! number of contour levels
  244. c     real*4 z(1:nc)             ! contour levels in increasing order
  245. c
  246.       subroutine conrec(d,ilb,iub,jlb,jub,x,y,nc,z)
  247.       real*4 d(ilb:iub,jlb:jub)
  248.       integer ilb,iub,jlb,jub
  249.       real*4 x(ilb:iub)
  250.       real*4 y(jlb:jub)
  251.       integer nc
  252.       real*4 z(1:nc)
  253. c
  254. c     Local declarations
  255. c
  256.       real*4 h(0:4)
  257.       integer sh(0:4)
  258.       real*4 xh(0:4),yh(0:4)
  259.       integer im(1:4),jm(1:4)
  260.       integer case
  261.       integer castab(-1:1,-1:1,-1:1)
  262.       integer p1,p2
  263. c
  264. c     Data
  265. c
  266.       data im/0,1,1,0/
  267.       data jm/0,0,1,1/
  268.       data castab/0,0,9,  0,1,5,  7,4,8,
  269.      1            0,3,6,  2,3,2,  6,3,0,
  270.      2            8,4,7,  5,1,0,  9,0,0/
  271. c
  272. c     Use statement functions for the line intersections
  273. c
  274.       xsect(p1,p2) = (h(p2)*xh(p1)-h(p1)*xh(p2))/(h(p2)-h(p1))
  275.       ysect(p1,p2) = (h(p2)*yh(p1)-h(p1)*yh(p2))/(h(p2)-h(p1))
  276. c
  277. c     Scan the arrays, top down, left to right within rows
  278. c
  279. 20    do 100 j=jub-1,jlb,-1
  280.          do 90 i=ilb,iub-1
  281.             dmin = min(d(i,j),d(i,j+1),d(i+1,j),d(i+1,j+1))
  282.             dmax = max(d(i,j),d(i,j+1),d(i+1,j),d(i+1,j+1))
  283.             if (dmax.ge.z(1) .and. dmin.le.z(nc)) then
  284.                do 80 k=1,nc
  285.                   if (z(k).ge.dmin .and. z(k).le.dmax) then
  286.                      do 22 m=4,0,-1
  287.                         if (m.gt.0) then
  288.                            h(m)=d(i+im(m),j+jm(m))-z(k)
  289.                            xh(m)=x(i+im(m))
  290.                            yh(m)=y(j+jm(m))
  291.                         else
  292.                            h(0)=0.25*(h(1)+h(2)+h(3)+h(4))
  293.                            xh(0)=0.5*(x(i)+x(i+1))
  294.                            yh(0)=0.5*(y(j)+y(j+1))
  295.                         endif
  296.                         if (h(m).gt.0.0) then
  297.                            sh(m)=+1
  298.                         else if (h(m).lt.0.0) then
  299.                            sh(m)=-1
  300.                         else
  301.                            sh(m)=0
  302.                         endif
  303. 22                   continue
  304. c
  305. c     Note: at this stage the relative heights of the corners and the
  306. c     centre are in the h array, and the corresponding coordinates are
  307. c     in the xh and yh arrays. The centre of the box is indexed by 0
  308. c     and the 4 corners by 1 to 4 as shown below.
  309. c     Each triangle is then indexed by the parameter m, and the 3
  310. c     vertices of each triangle are indexed by parameters m1,m2,and m3.
  311. c     It is assumed that the centre of the box is always vertex 2 though
  312. c     this isimportant only when all 3 vertices lie exactly on the same
  313. c     contour level, in which case only the side of the box is drawn.
  314. c
  315. c
  316. c           vertex 4 +-------------------+ vertex 3
  317. c                    | \               / |
  318. c                    |   \    m-3    /   |
  319. c                    |     \       /     |
  320. c                    |       \   /       |
  321. c                    |  m=2    X   m=2   |       the centre is vertex 0
  322. c                    |       /   \       |
  323. c                    |     /       \     |
  324. c                    |   /    m=1    \   |
  325. c                    | /               \ |
  326. c           vertex 1 +-------------------+ vertex 2
  327. c
  328. c
  329. c
  330. c                    Scan each triangle in the box
  331. c
  332.                      do 60 m=1,4
  333.                         m1=m
  334.                         m2=0
  335.                         if (m.ne.4) then
  336.                            m3=m+1
  337.                         else
  338.                            m3=1
  339.                         endif
  340.                         case = castab(sh(m1),sh(m2),sh(m3))
  341.                         if (case.ne.0) then
  342.                            goto (31,32,33,34,35,36,37,38,39),case
  343. c
  344. c     Case 1 - Line between vertices 1 and 2
  345. c
  346. 31                            x1=xh(m1)
  347.                               y1=yh(m1)
  348.                               x2=xh(m2)
  349.                               y2=yh(m2)
  350.                               goto 40
  351. c
  352. c     Case 2 - Line between vertices 2 and 3
  353. c
  354. 32                            x1=xh(m2)
  355.                               y1=yh(m2)
  356.                               x2=xh(m3)
  357.                               y2=yh(m3)
  358.                               goto 40
  359. c
  360. c     Case 3 - Line between vertices 3 and 1
  361. c
  362. 33                            x1=xh(m3)
  363.                               y1=yh(m3)
  364.                               x2=xh(m1)
  365.                               y2=yh(m1)
  366.                               goto 40
  367. c
  368. c     Case 4 - Line between vertex 1 and side 2-3
  369. c
  370. 34                            x1=xh(m1)
  371.                               y1=yh(m1)
  372.                               x2=xsect(m2,m3)
  373.                               y2=ysect(m2,m3)
  374.                               goto 40
  375. c
  376. c     Case 5 - Line between vertex 2 and side 3-1
  377. c
  378. 35                            x1=xh(m2)
  379.                               y1=yh(m2)
  380.                               x2=xsect(m3,m1)
  381.                               y2=ysect(m3,m1)
  382.                               goto 40
  383. c
  384. c     Case 6 - Line between vertex 3 and side 1-2
  385. c
  386. 36                            x1=xh(m3)
  387.                               y1=yh(m3)
  388.                               x2=xsect(m1,m2)
  389.                               y2=ysect(m1,m2)
  390.                               goto 40
  391. c
  392. c     Case 7 - Line between sides 1-2 and 2-3
  393. c
  394. 37                            x1=xsect(m1,m2)
  395.                               y1=ysect(m1,m2)
  396.                               x2=xsect(m2,m3)
  397.                               y2=ysect(m2,m3)
  398.                               goto 40
  399. c
  400. c     Case 8 - Line between sides 2-3 and 3-1
  401. c
  402. 38                            x1=xsect(m2,m3)
  403.                               y1=ysect(m2,m3)
  404.                               x2=xsect(m3,m1)
  405.                               y2=ysect(m3,m1)
  406.                               goto 40
  407. c
  408. c     Case 9 - Line between sides 3-1 and 1-2
  409. c
  410. 39                            x1=xsect(m3,m1)
  411.                               y1=ysect(m3,m1)
  412.                               x2=xsect(m1,m2)
  413.                               y2=ysect(m1,m2)
  414.                               goto 40
  415. 40                         call vecout(x1,y1,x2,y2,z(k))
  416.                         endif
  417. 60                   continue
  418.                   endif
  419. 80             continue
  420.             endif
  421. 90       continue
  422. 100   continue
  423.       return
  424.       end
  425. c
  426. c======================================================================
  427. c
  428. c     This is a sample vector output routine. For a local environment
  429. c     either replace the VECOUT call in the main line, or better
  430. c     replace the contents of this subroutine between the *'s shown.
  431. c
  432. c     There is often the requirement to distinguish each contour
  433. c     line with a different colour or a different line style. This
  434. c     can be done in many ways using the contour values z for a
  435. c     particular line segment.
  436. c
  437.       subroutine vecout(x1,y1,x2,y2,z)
  438.       implicit none
  439.       real*4 x1,y1,x2,y2,z
  440. c
  441. c***** Replace from here
  442. c
  443. c     The following should be ignored since it is specific to
  444. c     the version of FORTRAN running on the Macintosh microcomputer
  445. c     on which this particular example was written.
  446. c
  447.       INTEGER LINETO
  448.       PARAMETER (LINETO=Z'89109000')
  449.       INTEGER MOVETO
  450.       PARAMETER (MOVETO=Z'89309000')
  451.       call toolbx(MOVETO,nint(x1),nint(y1))
  452.       call toolbx(LINETO,nint(x2),nint(y2))
  453. c
  454. c***** To here
  455. c
  456.       return
  457.       end
  458.  
  459. ****** IBM BASIC CODE *******
  460.  
  461. 1000 REM 
  462. ====================================================================
  463. 1010 REM Documentation
  464. 1020 REM =============
  465. 1030 REM
  466. 1040 REM Originally written in FORTRAN-77
  467. 1050 REM Translated to BASICA on the IBM-PC microcomputer
  468. 1060 REM By Paul D. Bourke, August 1987
  469. 1070 REM
  470. 1080 REM Input variables to CONREC
  471. 1090 REM    d(0:iub,0:jub)   matrix for the data surface height values
  472. 1100 REM    iub, jub         index bounds of the data array
  473. 1110 REM    x(0:iub)         data array for column coordinates
  474. 1120 REM    y(0:jub)         data array for row coordinates
  475. 1130 REM    nc               number of contour levels
  476. 1140 REM    z(0:nc-1)        contour levels in increasing order
  477. 1150 REM    It is assumed that the logical variables TRUE and FALSE are  
  478. 1160 REM    either available or have been defined in the main line.
  479. 1170 REM
  480. 1180 REM NOTE: All arrays have a minimum index of 0
  481. 1190 REM
  482. 1200 REM Local declarations for CONREC
  483. 1210 REM =============================
  484. 1220 DIM H(4)     '    relative heights of the box above contour
  485. 1230 DIM ISH(4)   '    sign of h()
  486. 1240 DIM XH(4)    '    x coordinates of box
  487. 1250 DIM YH(4)    '    y coordinates of box
  488. 1260 DIM IM(3)    '    mapping from vertex numbers to x offsets
  489. 1270 IM(0)=0 : IM(1)=1 : IM(2)=1 : IM(3)=0
  490. 1280 DIM JM(3)    '    mapping from vertex numbers to y offsets
  491. 1290 JM(0)=0 : JM(1)=0 : JM(2)=1 : JM(3)=1
  492. 1300 DIM CASTAB(2,2,2)
  493. 1310 DATA 0,0,8,0,2,5,7,6,9, 0,3,4,1,3,1,4,3,0, 9,6,7,5,2,0,8,0,0
  494. 1320 FOR K=0 TO 2
  495. 1330    FOR J=0 TO 2
  496. 1340       FOR I=0 TO 2
  497. 1350          READ CASTAB(K,J,I)
  498. 1360       NEXT I
  499. 1370    NEXT J
  500. 1380 NEXT K
  501. 1390 REM
  502. 1400 REM Check the input parameters for validity
  503. 1410 REM PRMERR should be tested by the calling routine and
  504. 1420 REM MSG$ printed if PRMERR is TRUE
  505. 1430 REM
  506. 1440 PRMERR = FALSE
  507. 1450 IF (IUB<=0 OR JUB<=0) THEN PRMERR = TRUE
  508. 1460 IF (NC<=0) THEN PRMERR = TRUE
  509. 1470 FOR K=1 TO NC-1
  510. 1480    IF (Z(K)<=Z(K-1)) THEN PRMERR = TRUE
  511. 1490 NEXT K
  512. 1500 IF (PRMERR) THEN MSG$ = "Error in input parameters" : RETURN
  513. 1510 REM
  514. 1520 REM Scan the array, top down, left to right
  515. 1530 REM =======================================
  516. 1540 FOR J=JUB-1 TO 0 STEP -1
  517. 1550 FOR I=0 TO IUB-1
  518. 1560 REM Find the lowest vertex in the rectangle
  519. 1570 IF (D(I,J)<D(I,J+1)) THEN DMIN = D(I,J) ELSE DMIN = D(I,J+1)
  520. 1580 IF (D(I+1,J)<DMIN)   THEN DMIN = D(I+1,J)
  521. 1590 IF (D(I+1,J+1)<DMIN) THEN DMIN = D(I+1,J+1)
  522. 1600 REM Find the highest vertex in the rectangle
  523. 1610 IF (D(I,J)>D(I,J+1)) THEN DMAX = D(I,J) ELSE DMAX = D(I,J+1)
  524. 1620 IF (D(I+1,J)>DMAX)   THEN DMAX = D(I+1,J)
  525. 1630 IF (D(I+1,J+1)>DMAX) THEN DMAX = D(I+1,J+1)
  526. 1640 IF (DMAX<Z(0) OR DMIN>Z(NC-1)) THEN GOTO 2700
  527. 1650 REM
  528. 1660 REM Draw each contour within this box
  529. 1670 FOR K=0 TO NC-1
  530. 1680 IF ((Z(K)<DMIN) OR (Z(K)>DMAX)) THEN GOTO 2690
  531. 1690 FOR M=4 TO 0 STEP -1
  532. 1700    IF M>0 THEN GOTO 1720
  533. 1710    IF M=0 THEN GOTO 1770
  534. 1720    REM m > 0
  535. 1730       H(M)  = D(I+IM(M-1),J+JM(M-1))-Z(K)
  536. 1740       XH(M) = X(I+IM(M-1))
  537. 1750       YH(M) = Y(J+JM(M-1))
  538. 1760       GOTO 1820
  539. 1770    REM m = 0
  540. 1780       H(0)  = (H(1)+H(2)+H(3)+H(4))/4
  541. 1790       XH(0) = (X(I)+X(I+1))/2
  542. 1800       YH(0) = (Y(J)+Y(J+1))/2
  543. 1810       GOTO 1820
  544. 1820    IF H(M) > 0 THEN ISH(M) = 2
  545. 1830    IF H(M) = 0 THEN ISH(M) = 1
  546. 1840    IF H(M) < 0 THEN ISH(M) = 0
  547. 1850 NEXT M
  548. 1860 REM
  549. 1870 REM Scan each triangle in the box
  550. 1880 FOR M=1 TO 4
  551. 1890    M1 = M
  552. 1900    M2 = 0
  553. 1910    M3 = M+1
  554. 1920    IF (M3 = 5) THEN M3 = 1
  555. 1930 REM
  556. 1940 REM Select the type of line/plane intersection
  557. 1950 CASE = CINT(CASTAB(ISH(M1),ISH(M2),ISH(M3)))
  558. 1960 IF CASE = 0 THEN GOTO 2680
  559. 1970 IF CASE = 1 THEN GOTO 2060
  560. 1980 IF CASE = 2 THEN GOTO 2120
  561. 1990 IF CASE = 3 THEN GOTO 2180
  562. 2000 IF CASE = 4 THEN GOTO 2240
  563. 2010 IF CASE = 5 THEN GOTO 2300
  564. 2020 IF CASE = 6 THEN GOTO 2360
  565. 2030 IF CASE = 7 THEN GOTO 2420
  566. 2040 IF CASE = 8 THEN GOTO 2480
  567. 2050 IF CASE = 9 THEN GOTO 2540
  568. 2060 REM Case 1: Line between vertices m1 and m2
  569. 2070    X1 = XH(M1)
  570. 2080    Y1 = YH(M1)
  571. 2090    X2 = XH(M2)
  572. 2100    Y2 = YH(M2)
  573. 2110    GOTO 2650
  574. 2120 REM Case 2: Line between vertices m2 and m3
  575. 2130    X1 = XH(M2)
  576. 2140    Y1 = YH(M2)
  577. 2150    X2 = XH(M3)
  578. 2160    Y2 = YH(M3)
  579. 2170    GOTO 2650
  580. 2180 REM Case 3: Line between vertices m3 and m1
  581. 2190    X1 = XH(M3)
  582. 2200    Y1 = YH(M3)
  583. 2210    X2 = XH(M1)
  584. 2220    Y2 = YH(M1)
  585. 2230    GOTO 2650
  586. 2240 REM Case 4: Line between vertex m1 and side m2-m3
  587. 2250    X1 = XH(M1)
  588. 2260    Y1 = YH(M1)
  589. 2270    X2 = (H(M3)*XH(M2)-H(M2)*XH(M3))/(H(M3)-H(M2))
  590. 2280    Y2 = (H(M3)*YH(M2)-H(M2)*YH(M3))/(H(M3)-H(M2))
  591. 2290    GOTO 2650
  592. 2300 REM Case 5: Line between vertex m2 and side m3-m1
  593. 2310    X1 = XH(M2)
  594. 2320    Y1 = YH(M2)
  595. 2330    X2 = (H(M1)*XH(M3)-H(M3)*XH(M1))/(H(M1)-H(M3))
  596. 2340    Y2 = (H(M1)*YH(M3)-H(M3)*YH(M1))/(H(M1)-H(M3))
  597. 2350    GOTO 2650
  598. 2360 REM Case 6: Line between vertex m3 and side m1-m2
  599. 2370    X1 = XH(M3)
  600. 2380    Y1 = YH(M3)
  601. 2390    X2 = (H(M2)*XH(M1)-H(M1)*XH(M2))/(H(M2)-H(M1))
  602. 2400    Y2 = (H(M2)*YH(M1)-H(M1)*YH(M2))/(H(M2)-H(M1))
  603. 2410    GOTO 2650
  604. 2420 REM CASE 7: Line between sides m1-m2 and m2-m3
  605. 2430    X1 = (H(M2)*XH(M1)-H(M1)*XH(M2))/(H(M2)-H(M1))
  606. 2440    Y1 = (H(M2)*YH(M1)-H(M1)*YH(M2))/(H(M2)-H(M1))
  607. 2450    X2 = (H(M3)*XH(M2)-H(M2)*XH(M3))/(H(M3)-H(M2))
  608. 2460    Y2 = (H(M3)*YH(M2)-H(M2)*YH(M3))/(H(M3)-H(M2))
  609. 2470    GOTO 2650
  610. 2480 REM Case 8: Line between sides m2-m3 and m3-m1
  611. 2490    X1 = (H(M3)*XH(M2)-H(M2)*XH(M3))/(H(M3)-H(M2))
  612. 2500    Y1 = (H(M3)*YH(M2)-H(M2)*YH(M3))/(H(M3)-H(M2))
  613. 2510    X2 = (H(M1)*XH(M3)-H(M3)*XH(M1))/(H(M1)-H(M3))
  614. 2520    Y2 = (H(M1)*YH(M3)-H(M3)*YH(M1))/(H(M1)-H(M3))
  615. 2530    GOTO 2650
  616. 2540 REM Case 9: Line between sides m3-m1 and m1-m2
  617. 2550    X1 = (H(M1)*XH(M3)-H(M3)*XH(M1))/(H(M1)-H(M3))
  618. 2560    Y1 = (H(M1)*YH(M3)-H(M3)*YH(M1))/(H(M1)-H(M3))
  619. 2570    X2 = (H(M2)*XH(M1)-H(M1)*XH(M2))/(H(M2)-H(M1))
  620. 2580    Y2 = (H(M2)*YH(M1)-H(M1)*YH(M2))/(H(M2)-H(M1))
  621. 2590 REM This is where a contour line segment of value z(i)
  622. 2600 REM is drawn between points (x1,y1) and (x2,y2)
  623. 2610 REM The exact details can vary depending on the particular
  624. 2620 REM graphics capabilities available. 
  625. 2630 REM If colour is available then each contour level could 
  626. 2640 REM be drawn in a new colour
  627. 2650 line (x1,y1)-(x2,y2)
  628. 2660 NEXT M
  629. 2670 NEXT K
  630. 2680 NEXT I
  631. 2690 NEXT J
  632. 2700 RETURN
  633. 2710 REM ===================================================================
  634.